!-------------------------------------------------------------------------------------
subroutine diagnostics( t, t_iter, Courant_in_z, &               ! inputs
                        OmegaBounce, SpeckleWidth,  &            ! inputs
                        dfu, f0vec, &                      ! inputs
                        Ezinternal, Ez2D, Ex2D, phiTot, phiInt)  ! inputs

!-------------------------------------------------------------------------------------

  use dimensions

  use fft_wrap, only   : fft_forward, fft_backward
  use parameters, only : kz, pi, iself, ny

  implicit none

  !-- global variables --

  real, intent(in)       :: t, t_iter, Courant_in_z, OmegaBounce, SpeckleWidth

  real, intent(in)       :: dfu(jz, jv, jx, jDirection)
  real, intent(in)       :: f0vec(jv)

  real, intent(in)       :: Ezinternal(jx,jz)
  real, intent(in)       :: Ez2D(jx,jz), Ex2D(jx,jz)

  complex, intent(in)    :: phiTot, phiInt

  !-- local variables --

  real       :: PhasePhi
  real, save :: PhasePhiOld
  real       :: KE, PotE, PotEx, potEz, Phi0Energy
  real       :: nTrapped, fsq


  complex    :: epsilon
  complex    :: EzFtInternal, EzFTtotal

            
  real       :: funiform(jv)

  complex    :: Phi1Dz(jz)


  integer    :: ix, iz, iv, iDirection  ! counters


  !-----------------------------------------------------------------------------------
 
      !phiTot = phi0T*kz*kz
      !phiInt = 2*density1DFT(1+ny)
                 
      if (iself==1) then
         epsilon = phiTot / (phiInt + phiTot)
         !phi = 2*density1DFT(1+ny)/kz**2 + phi0T
       else
         epsilon = ( phiTot - phiInt ) / phiTot 
         ! = phiExternal/phiTotal = (phiTotal-phiInternal)/phiTotal
      end if
      
      nTrapped=0.0
      fsq=0
        
        
      do iv=1,jv
         funiform(iv)=0.0
         do ix=1,jx
            do iz=1,jz
               do iDirection = 1,2
                  funiform(iv) = funiform(iv)+dfu(iz,iv,ix,iDirection)/(jx*jz)
                  fsq = fsq + ( f0vec(iv) + dfu(iz,iv,ix,iDirection) )**2 * dx*dz*dv
               enddo! iDir
            enddo ! z
         enddo !  x
      end do !  iv
        
        
      nTrapped = nTrapped/jz
        
      PotE=0;  PotEx=0; potEz=0; PHi0Energy=0
      do iz=1,jz
         do ix=1,jx
            PotEx = PotEx + 0.5 * Ex2D(ix,iz)**2 /(jx*jz)
            potEz = potEz + 0.5 * Ez2D(ix,iz)**2 /(jx*jz)
         enddo
      enddo
        
      potE = PotEx + potEz
        
      KE=0

      do iv=1,jv 
         KE =  KE + funiform(iv) * 0.5 * vvec(iv)**2 * dv
      enddo
        
      Phi1Dz(:) = Ezinternal(jx/2 + 1, :)
      call fft_forward(jz,Phi1Dz,Phi1Dz)
      Phi1Dz=Phi1Dz/float(jz)
      EzFtInternal = Phi1Dz(1+ny)
        
      Phi1Dz(:) = Ez2D(jx/2 + 1, :)
      call fft_forward(jz,Phi1Dz,Phi1Dz)
      Phi1Dz=Phi1Dz/float(jz)
      EzFTtotal = Phi1Dz(1+ny)

      PhasePhi = -aimag(log(EzFtInternal)) ! really phase of E
        
      do while ( ( PhasePhi - PhasePhiOld ) <  - pi )
         PhasePhi = PhasePhi + 2*pi
      enddo
      do while ( ( PhasePhi - PhasePhiOld ) >  + pi )
         PhasePhi = PhasePhi - 2*pi
      enddo
        
      PhasePhiOld = PhasePhi   

     !------------------------------
       

     ! t_iter =  float(it-itMin)/nt

      open(3, file='vstime.txt', position='append')
        
      if (iself==1) then
         write(3,501) t, 2*abs(EzFtInternal) / kz, 2*abs(EzFTtotal) / kz, PhasePhi, KE, potEx, potEz, &
              OmegaBounce*SpeckleWidth*0.37233, Courant_in_z
         write(*,501) t_iter, t, 2*abs(EzFtInternal) / kz, &
              OmegaBounce*SpeckleWidth*0.37233, potEx, potEz, Courant_in_z
      else
         write(3,501) t, 1/cabs(epsilon), real(epsilon), aimag(epsilon), &
              kz*kz * real( 1 - epsilon ), &
              kz*kz * aimag( 1 - epsilon ), &
              KE,  PHi0Energy, KE+PHi0Energy, fsq, nTrapped
         write(*,501) t_iter, kz*t, real(epsilon), KE+PHi0Energy, KE
      endif

      close(3)

    501   format(  6400(e12.6,2x)  )


end subroutine diagnostics

!-------------------------------------------------------------------------------------
subroutine write_csi()
!-------------------------------------------------------------------------------------

  use parameters, only : phi0, kz, vPhase

  implicit none

  complex    :: epsilon, epsilon0, Csi0, dCsi0dv, dEpsilonDomega
 

      !------------------------------

      open(5, file='output', position='append')
                                                                             
      call getCsi0(Csi0,cmplx(vPhase,0.0),dCsi0dv)

      write(5,*) 'Re(Csi0)= ', real(Csi0), &
           '    Imag(Csi0/kz^2)=', aimag(Csi0)/kz**2, &
           '  dCsi0dv=',dCsi0dv

      write(5,*); 
      write(5,*)' dCsi0dv =', dCsi0dv
      if(  real(Csi0)  >  0  ) then
         epsilon0 = 1 - Csi0/ real(Csi0)
         dEpsilonDomega = - dCsi0dv / sqrt(real(Csi0))**3
         write(5,*) ' kz at resonance=', sqrt(real(Csi0)), &
              '   Landau damping=', aimag( epsilon0 ) / real( dEpsilonDomega  )
         write(5,*)'   Landau damping=', aimag( epsilon0 ) / real( dEpsilonDomega  )
      endif


      call getCsi0(Csi0,cmplx(vPhase,0.0),dCsi0dv)
      write(5,*) 'Re(Csi0)= ', real(Csi0), &
           '    Imag(Csi0/kz^2)=', aimag(Csi0)/kz**2, &
           '  dCsi0dv=',dCsi0dv
      write(5,*); 
      write(5,*)' dCsi0dv =', dCsi0dv
      epsilon0 = 1 - Csi0/ kz**2

      write(5,*) '   |PhiInternal|=', phi0*abs(   1/epsilon0 - 1 )

      close(5)


end subroutine write_csi

!-------------------------------------------------------------------------------------
subroutine write_fdist(dfu, OmegaBounce)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only : kz, vPhase, iWaveFrame, pi

  implicit none


  real, intent(in)       :: dfu(jz, jv, jx, jDirection)
  real, intent(in)       :: OmegaBounce  

 
  integer    :: ix, iz, iv
  real       :: v, f0


     f0(v)=( 1/sqrt(2*pi) )*exp(-.5*(v+vPhase*iWaveFrame)*(v+vPhase*iWaveFrame))

     open(16, file='df_vs_zv.dat')
     open(18, file='f_vs_zv.dat')
     open(13,file='f_trapped_vs_zv.dat')

     write(16,*) '% v\z';  write(16,*)  'NA', (zvec(iz),iz=1,jz)
     write(18,*) '% v\z';  write(18,*)  'NA', (zvec(iz),iz=1,jz)
     write(13,*) '% v\z';  write(13,*)  'NA', (zvec(iz),iz=1,jz)

     ix=jx/2 + 1    
     do iv=1,jv
        write(16,501) vvec(iv),  (    dfu(iz,iv,ix,1) + dfu(iz,iv,ix,2)               , iz=1,jz   ) ! dfzv
        write(18,501) vvec(iv),  (    dfu(iz,iv,ix,1) + dfu(iz,iv,ix,2) + f0(vvec(iv)), iz=1,jz   ) ! fzv
        if  (   abs( vvec(iv) ) < sqrt (2 * OmegaBounce/kz )    )then                  
           write(13,501) vvec(iv), ( dfu(iz,iv,ix,1) + dfu(iz,iv,ix,2) + f0(vvec(iv)), iz=1,jz   ) ! fzv trapped
        endif
     end do;

     close(13)
     close(18)
     close(16)

 501   format(  6400(e12.6,2x)  )


end subroutine write_fdist


!-------------------------------------------------------------------------------------
subroutine write_means(dfu, fvc0T)
!-------------------------------------------------------------------------------------

  use dimensions
  use fft_wrap, only   : fft_forward, fft_backward


  implicit none


  real, intent(in)       :: dfu(jz, jv, jx, jDirection)
  complex, intent(in)    :: fvc0T(jv)
 
  real       :: MeanFTvSq(jv), MeanFTkSq(jz)
  complex    :: fz(jz)
  complex    :: fvc(jv), fvc0(jv)

  integer    :: ix, iz, iv, ikz, ip
  real       :: MySum




      MeanFTvSq(:)=0; MeanFTkSq(:)=0

      ix=jx/2
      do iv=1,jv         
         fz(:) = dfu(:,iv,ix,1) + dfu(:,iv,ix,2)
         call fft_forward(jz,fz,fz)  
         fz=fz/float(jz)
         call dealias(nkz,fz)
         MeanFTkSq(:)=MeanFTkSq(:)+abs(fz(:)) **2                        
      end do
 
      ix=jx/2 
      do iz=1,jz 
         fvc(:) = dfu(iz,:,ix,1) + dfu(iz,:,ix,2)       
         call fft_forward(jv,fvc,fvc)
         fvc=fvc/float(jv)
         call dealias(jv,fvc)
         MeanFTvSq=MeanFTvSq+abs(fvc) **2
      end do

          
      open(55,file='FTvSQ_avg.txt')
      write(55,*)'% p   <FTvSQ>   F0FT^2   Cum(p^4*<FTvSQ>)'; MySum=0
      do ip=1,jv/2  - jv/6 -1
         MySum = MySum + pvec(ip)**4 * MeanFTvSq(ip)
         write(55,501)pvec(ip),MeanFTvSq(ip), abs(fvc0T(ip))**2, MySum
      end do
      close(55)
      
      open(55,file='FTkSQ_avg.txt')
      write(55,*)'% kz   <FTkSQ>   Cum(kz^4*<FTkSQ>) ' ;  MySum=0
      do ikz=1,nkz/2  - nkz/6 -1
         MySum=MySum + kZvec(ikz)**4 * MeanFTkSq(ikz)
         write(55,501)kZvec(ikz),MeanFTkSq(ikz)    , MySum
      end do
      close(55)
       


    501   format(  6400(e12.6,2x)  )


end subroutine write_means


!-------------------------------------------------------------------------------------
subroutine write_umat_spectral(uzMat, uxMat)              ! inputs 
!-------------------------------------------------------------------------------------

  use dimensions
  use fft_wrap, only   : fft_forward, fft_backward
  use parameters, only : u

  implicit none

  real, intent(in)        :: uzMat(jx, jz, jDirection)
  real, intent(in)        :: uxMat(jx, jz, jDirection)

  integer    :: ix, iz, iv, ikz, ikx, iDirection


  real       :: avUx(jx)
  complex    :: du(jx)
  complex    :: dummy(jx)


  avUx = 0 !  real avUx(jx)    DIAGNOSTIC
  iDirection = 1
  do iz=1,jz
     du(:) = uxMat(:,iz,iDirection)
     call fft_forward(jx,du,du)
     du=du/float(jx)
     avUx = avUx + abs(du)**2 / jz
  enddo  ! iz

  open(12,file="Ux_kx_Sq.txt")
  write(12,*)"%kx  <|Ux_kx_Sq|^2>"      
  do ikx=1,nkx/2  - nkx/6 -1
     write(12,501) kxVec(ikx), avUx(ikx)
  end do
  close(12)


  avUx = 0 !  real avUx(jx)    DIAGNOSTIC
  iDirection = 1
  do iz=1,jz
     du(:) = uzMat(:,iz,iDirection)
     call fft_forward(jx,du,du)
     du=du/float(jx)
     avUx = avUx + abs(du)**2 / jz
  enddo


  open(12,file="Uz_kx_Sq.txt")
  write(12,*)"%kx  <|Uz_kx_Sq|^2>"      
  do ikx=1,nkx/2  - nkx/6 -1
     write(12,501) kxVec(ikx), avUx(ikx)
  end do; 
  close(12)



  iz = jz/3
  dummy(:) = uxMat(:,iz,1) - u
  call fft_forward(jx,dummy,dummy)

  open(26,file='uxPlus_vs_kx.txt')
  write(26,*)'% kx   |Uplus(kx,iz=jz/3)|^2'
  do ix = 1, jx/2
     write(26,501) kXvec(ix), abs(dummy(ix))**2 
  enddo
  close(26)

           

501  format(  6400(e12.6,2x)  )


end subroutine write_umat_spectral

    
!-------------------------------------------------------------------------------------
! subroutine getCsi0
!-------------------------------------------------------------------------------------

  subroutine getCsi0(Csi0,vc,dCsi0dv)     !	real(vc)>=0,  aimag(vc)>=0, or both <0
  implicit none
  complex Csi0,vc,dCsi0dv         !	dCsi/dzeta
  real x,mu
  complex Z,Zprime,Zprime2
  
  x=abs ( real(vc) )
  mu=abs( aimag(vc) )
  
  call  getZstuff(x/1.41421,Z,Zprime,Zprime2)     !	1.41421=sqrt(2)
  
  dCsi0dv=3.53553E-01*Zprime2     !	3.53553E-01=1/2**1.5
  
  Csi0=0.5*Zprime + (0.0,1.0)*mu*dCsi0dv
  
  return
  end 


!-------------------------------------------------------------------------------------
! subroutine getZstuff
!-------------------------------------------------------------------------------------

  subroutine getZstuff(x,Z,Zprime,Zprime2)        ! should be  good  for  x pos or neg
  implicit none
  real x,pi, DAWS, MyDAWS !  DAWS is Dawson integral, as give by IMSL library.
  complex Z,i,Zprime,Zprime2
  data pi,i/3.14159265,(0.0,1.0)/
  
!  Z=i*sqrt(pi)*exp(-x*x) -2*DAWS(x)
   Z=i*sqrt(pi)*exp(-x*x) -2*MyDaws(x)  !  use this if DAWS is unavailable
  
  Zprime=-2*(1+x*Z) !		dZ/dz
  Zprime2=-2*x*Zprime-2*Z     !	d2Z/dx2
  return
  end
  
  real function MyDaws(z1)  !  reproduces DAWS to high accuracy
  implicit none
  real x,z1,Z
  x=abs(z1)
   if( (0.0.le.x).and.(x.lt.1.0) )Z=-0.6953*x**4 + 1.692*x**3  -0.07885*x**2  -1.994*x
   if( (1.0.le.x).and.(x.lt.2.0) )Z=+0.2535577*x**4  -1.905759*x**3 + 5.158011*x**2  -5.463602*x+ 0.8816458
   if( (2.0.le.x).and.(x.lt.2.5) )Z=-(1/x)*( 1+0.5/x**2 +0.75/x**4 +(15./8.)/x**6 ) &
&       +1.563421E-1*x**4  -1.516732*x**3 + 5.520618*x**2  -8.926533*x + 5.400882
   if( x.ge.2.5 )Z=-(1/x)*( 1+0.5/x**2 +0.75/x**4 +(15./8.)/x**6 )  -6.124581*Exp(-3.024965*x)

   MyDaws=-0.5*Z
   if(  z1 < 0  )MyDaws=-MyDaws
   return
   end 
  
  

!-------------------------------------------------------------------------------------
! subroutine write2D_vs_xz
!-------------------------------------------------------------------------------------

  subroutine write2D_vs_xz(fname, data)

    use dimensions

    implicit none

    real          :: data(jx,jz)
    character (*) :: fname
    

    integer ix, iz

    open(206,file=fname);

    write(206,*)'% z\x'; 
    write(206,*) 'NA', (  xvec(ix), ix=1,jx )
    do iz = 1, jz
      write(206,501) zvec(iz), ( data(ix,iz), ix=1,jx )
    enddo

    close(206)

    501   format(  6400(e12.6,2x)  )

  end

!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------


